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PACS 64 . 70 . Fx - Liquid- vapor transitions 

PACS 47 . 1 1 . Qr - Lattice gas 

PACS 64 . 70 . Md - Transitions in liquid crystals 

Abstract. - We present Monte Carlo simulation results of the two-dimensional Zwanzig fluid, 
which consists of hard line segments which may orient either horizontally or vertically. At a 
certain critical fugacity, we observe a phase transition with a two-dimensional Ising critical point. 
Above the transition point, the system is in an ordered state, with the majority of particles being 
either horizontally or vertically aligned. In contrast to previous work, we identify the transition 
as being of the liquid-gas type, as opposed to isotropic-to-nematic. This interpretation naturally 
accounts for the observed Ising critical behavior. Furthermore, when the Zwanzig fluid is extended 
to more allowed particle orientations, we argue that in some cases the symmetry of a q-state Potts 
model with q > 2 arises. This observation is used to interpret a number of previous results. 
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Introduction. — In a seminal paper [1], Onsager 
demonstrated that infinitely slender rods in three dimen- 
sions undergo a first-order isotropic-to-nematic (IN) tran- 
sition. In the nematic phase, there is long-ranged align- 
ment of the particles, while in the isotropic phase the par- 
ticle orientations are essentially random. In contrast, in 
two dimensions, long-ranged nematic order is generally ab- 
sent. For a certain class of liquid crystal pair potentials, 
the absence of nematic order can be proved rigorously [2] , 
while simulations using different potentials also indicate 
its absence in the thermodynamic limit [3,4]. Of course, 
these results do not imply that there can be no phase tran- 
sition in two-dimensional (2D) liquid crystals, but rather 
that any such transition does not lead to nematic order. 

Interestingly, a number of papers have appeared re- 
cently [5-8] in which the IN transition was studied in two 
dimensions. The transition was shown to belong to the 
universality class of the 2D Ising model. In accordance 
with the Ising model, this implies the formation of a finite 
nematic order parameter above the critical density, which 
seems to contrast the results of [3,4], where nematic or- 
der was found to vanish in the thermodynamic limit. The 
results of [5-8] thus raise a number of questions. First of 
all, how can we understand the formation of finite nematic 
order in these 2D systems and, secondly, what is the origin 
of the Ising critical point? In this work, these questions 



will be answered. 

The absence of nematic order in many two dimen- 
sional systems is a consequence of the Mermin- Wagner 
theorem [9-11]. As is well known, this theorem applies 
when the particle orientations are continuous. However, 
when the orientations become discretized, Mermin- Wagner 
no longer applies, and the corresponding phase behav- 
ior changes dramatically. A famous example of a liquid 
crystal model with discrete orientations is the Zwanzig 
model [12], where the particles are treated as rigid rods. 
The particle positions are continuous, but the molecular 
axis may only point in mutually perpendicular directions. 
In two dimensions, this implies a system of line segments, 
which may either point horizontally or vertically. The in- 
teractions are of the excluded volume type, meaning that 
particles may not overlap with each other. A further ap- 
proximation is to also make the particle positions discrete, 
i.e. to restrict the line segments to the sites of a square lat- 
tice, and to let each segment occupy k consecutive sites. 
This is precisely the model studied in [5-8]. For k = 2 
one recovers the dimer model [13], while for k — > oo one 
approaches the 2D Zwanzig model. Provided k > 7 [8], 
at some threshold density, one finds a transition to a "ne- 
matic" phase, in which most particles point either hori- 
zontally ( "A particles" ) or vertically ( "B particles" ) . 

However, is this transition truly an IN transition? In 
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this model, there is symmetry under the exchange of A 
and B particles. That is, given a valid configuration, 
i.e. one without overlaps, a new valid configuration can 
be obtained by replacing each A particle with a B particle 
and vice-versa. Under this operation, the order param- 
eter \Na — Nb\/(Na + Nb) [5] remains invariant, with 
Ni the number of particles of type i. Clearly, a liquid 
crystal with continuous orientations cannot exhibit this 
symmetry. The observed symmetry rather resembles the 
particle-hole symmetry of the lattice gas, or the up-down 
symmetry in the Ising model. This suggests that the IN 
transition of [5-8] is really an unmixing or liquid-gas tran- 
sition. This also accounts for the observed Ising critical 
behavior, since it is well known that liquid-gas and un- 
mixing transitions in fluids with short ranged interactions 
belong to this class. 

The interpretation in terms of a liquid-gas transition is 
also consistent with the original paper by Zwanzig [12]. 
Here, it was already mentioned that, in two dimensions, 
a mapping of the Zwanzig model onto a mixture of green 
and red squares [14] is possible, whereby squares of dif- 
ferent color may not overlap. This is, of course, just a 
variant of the Widom-Rowlinson mixture of spherical par- 
ticles in two dimensions [15], in which the existence of a 
liquid-gas transition is not debated. Note that the ori- 
gin of liquid-gas transitions in these systems stems from 
depletion. One could envision formally integrating out, 
say, the "red" species, yielding a one component fluid of 
"green" species, interacting via an effective short-ranged 
potential with attractive part [16]. Hence, these systems 
resemble simple fluids, such as the Lcnard- Jones fluid, and 
are expected to yield similar phase diagrams as a result. 

In this Letter, these ideas will be applied to the 2D 
Zwanzig model using computer simulation. We first spec- 
ify the 2D Zwanzig model, and describe the simulation 
method. Next, we show that the transition in this model 
indeed corresponds to a liquid-gas transition, rather than 
IN. In particular, we demonstrate that a binodal can be 
constructed, which terminates at an Ising critical point. 
We also provide estimates for the line-tension between co- 
existing domains in the two-phase (ordered) region of the 
phase diagram. Finally, we present a summary and de- 
tailed conclusion, where we emphasize that care must be 
taken when modeling liquid crystal phase transitions using 
only a discrete set of orientations. 

model and simulation method. We consider the 
2D Zwanzig model, which consists of infinitesimally thin 
hard rods of unit length. Since the interactions are hard- 
core, temperature does not play a role, and factors of ksT 
are set to unity throughout (with ks the Boltzmann con- 
stant and T the temperature). The rods may be aligned 
horizontally ("A particles") or vertically ("B particles"). 
The particle positions are confined to a periodic 2D square 
of area V. Since the rods are infinitesimally thin, there is 
only a hard-core interaction between A and B particles 
(which may thus not overlap). Hence, each A particle is 



surrounded by a depletion zone, which may not contain 
the centers of any B particles (note that the depletion 
zone for this model is just the unit square). We consider 
a grand canonical simulation ensemble, with the respec- 
tive chemical potentials fi and fig, of A and B particles, 
being the relevant thermodynamic parameters; the actual 
number of particles in the system is a fluctuating quantity. 
The aim of the simulations is to measure the distribution 
P(N\fx, /is), defined as the probability to observe a system 
containing TV particles of type A, at chemical potentials (i 
and fiB- 

During the simulations, insertion and removal of parti- 
cles are performed using a grand canonical cluster move 
[17,18]. With equal probability, we attempt to insert an A 
particle, or we attempt to remove one. When inserting, a 
single A particle is tentatively placed at a random location 
in the simulation box. This will generally lead to overlap 
with some B particles, say ub of them. The overlapping 
B particles are removed from the box, and the resulting 
state is accepted with probability 
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n B > A 
otherwise, 



with A a parameter to be specified later, and Nb the 
total number of B particles in the system at the beginning 
of the move. In the above, we have also introduced the 
respective fugacities z = cxp(/x) and Zb = exp(/xs), of A 
and B particles. During removal, one A particle is picked 
randomly and deleted, and n# centers of B particles are 
distributed randomly into the depletion zone of the just 
deleted A particle, with ub a uniform random number < 
ub < A. If any of the inserted B particles overlap with 
A particles, the move is rejected, otherwise it is accepted 
with probability 
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The reader may verify that this algorithm fulfills detailed 
balance [17, 18]. The factorials count the number of ways 
in which ns particles can be distributed onto the unit 
square. Note that the thermal wavelength has been set to 
unity for clarity. The parameter A must be high enough 
such that the insertion of A particles into a pure phase 
of B particles is efficient. For the present model, a pure 
phase of B particles is just an ideal gas, for which density 
equals fugacity. Hence, the depletion zone contains zb B 
particles on average, with Poissonian fluctuations. Conse- 
quently, A should somewhat exceed this value; we found 
that A = z B + ^/z~b + 2 gave good results. 

The phase transition in our model is characterized by 
a free energy functional featuring two minima separated 
by a barrier. Apart from a minus sign, the free energy 
is just the logarithm of the distribution P(N\fi, ub) that 
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we wish to find, and we define W(N) = hxP(N). Com- 
puter simulations which directly sample the Boltzmann 
distribution are not efficient then, since these tend to "get 
stuck" in one of the minima, and rarely cross the barrier. 
To overcome this problem, we combine the grand canon- 
ical cluster move with a biased sampling method called 
successive umbrella sampling (SUS) [19]. In SUS, W(N) 
is obtained recursively by splitting the simulation into a 
number of windows. In the first window, the number of 
A particles is allowed to fluctuate between and 1, in the 
second window between 1 and 2, and so forth. There is 
no restriction on the number of B particles though, and 
Ng fluctuates freely in each window. By simulating the 
7V-th window, one immediately obtains the free energy 
difference 

AF(N\fj,,fj, B ) = W(N) - W(N-1) = In (C+/C_) , (1) 

with C+ (C_) the number of times that the (unbiased) 
simulation was in a state with N (N—l) particles of type A 
(irrespective of iVs). Obviously, the free energy difference 
depends on \i and lib- Once the free energy differences 
have been measured over a range of windows, W(N) can 
be constructed via recursion 

W(0) = 0, W(N) = W(N - 1) + AF(N). (2) 

Using that P{N) oc e w ( N \ one trivially converts to the 
sought-for distribution P(N\li, llb)- 

An additional ingredient of this work is histogram 
reweighting, which we use to extrapolate simulation 
data obtained at {li,llb) to different chemical potentials 
(//, li' b )- In order to extrapolate in /is, we also require the 
distributions R(Nb\N, lib), defined as the probability to 
observe a system containing Nb particles of type B, when 
the number of A particles equals N, at chemical poten- 
tial lib (since R(Nb\N, lib) is obtained for fixed N, there 
is no dependence on ll). For example, R(Nb\0, lib) is the 
distribution in B particles when no A particles are present 
(this corresponds to an ideal gas at chemical potential li b , 
and hence a single Poissonian peak) . The expression to ex- 
trapolate the free energy difference obtained at (jx, lib) to 
(//, li' b ) then becomes 



AF(N\li',li' b ) = 
AF{N\li,li b ) 



In- 



Z(N) 



Z(N -I) 



with 



Z(N) 



En b R(Nb\N^ b ) 



(3) 



(4) 



Note that Z{N) is simply the relative change in the "vol- 
ume" of R(Nb\N, (j,b) when extrapolating from lib — * h'b- 
By using Eq.(3), in combination with the recursion rela- 
tion, it becomes possible to construct P{N\fx' , li' b ), with- 
out actually having to perform a simulation at /i' and \x' B . 



The quality of Eq.(3) deteriorates when the range in 
chemical potential over which one extrapolates becomes 
large. For each system size, we therefore perform a series 
of i = 1, . . . , k SUS simulations, over a range of chemical 
potentials [n and llbs- Due to symmetry, we set Li{ = 
/JiB,i for convenience, although this is not essential. Using 
Eq.(3), each one of these simulations yields an estimate of 
the free energy difference AFi(N\Li, /is) ± ffi, where Cj is 
an estimate of the statistical error. Next, we weight each 
estimate with its inverse square error to obtain the best 
estimate of the free energy difference 



AF boBt (N\n, hb) = 



(5) 



which is then fed into the recursion relation to construct 
the best estimate of P{N\ll 1 (j.b) for some [i, lib ofinterest. 
To derive o^, we note that statistical errors occur in the 
counts C+ and C_, as well as in the histogram entries 
R(Nb\N, hb)- In principle, these errors are Poissonian 

*[C+] = y/C^, a[C-] = y/CZ, 

<j[R(N b \N 7 lib)} = x/R(Nb\N,lib), 

but only if the data are normalized to the number of inde- 
pendent measurements. This requires knowledge of the 
correlation time r, which is computationally expensive 
to obtain, since the acceptance rate of the cluster move, 
and hence r, depend sensitively on density, composition, 
and system size. Instead, we follow a more pragmatic ap- 
proach, whereby C+ and C_ are normalized to the num- 
ber of accepted cluster moves in the window. Similarly, 
the histogram R(Nb\N, lib) is normalized to the number 
of cluster moves which resulted in a state with N particles 
of type A and which involved a change in the number of B 
particles. With these choices, we assume that the statisti- 
cal errors become Poissonian, i.e. proportional to square- 
roots, with a common proportionality constant. Next, a 
propagation of errors calculation can be applied to Eq.(3) 
and Eq.(4) to derive cr^. A final important optimization is 
to combine the set of histograms R{(Nb\N, llb,i) from the 
various SUS simulations into one best estimate using the 
multiple histogram method [20,21], and to subsequently 
use this best estimate to calculate Z{N) of Eq.(4). 

Results. — In Fig. 1, we show distributions W(N) = 
\nP(N), for a number of fugacities zb- For sufficiently 
large Zb, the distributions develop two pronounced peaks: 
one at low density p = N/V of A particles, and one at high 
density. Note that the bimodal structure only shows-up 
if the fugacity z of the A particles is chosen suitably. For 
the present model, of course, we may set z = zb due to 
symmetry, which was adopted throughout this work. For 
asymmetric fluids, choosing z is less straightforward, and 
many criteria can, in fact, be defined [22]. 

The connection to the liquid-gas transition becomes 
clear if one "identifies" the low-density peak with the gas 
phase, the high-density peak with the liquid, and zb with 
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Fig. 1: Distributions W(N) obtained in a 40 x 40 system. The 
distributions were measured at zb = 5.0,5.1,5.2,5.3,5.4, and 
clearly illustrate the formation of the double-peaked structure 
with increasing zb- 



inverse temperature. The region between the peaks re- 
flects phase coexistence, whereby both phases appear si- 
multaneously. In simulations, the coexistence can be visu- 
alized directly (Fig. 2). Note that a rectangular simulation 
box is used, and so the interfaces form parallel to the short 
edge, since this minimizes the total amount of interface 
(due to periodic boundary conditions, two interfaces are 
actually present). The corresponding distribution W(N) 
for the rectangular system is shown in Fig. 3. Note the flat 
region between the peaks, implying that interactions be- 
tween the interfaces are absent. Following Binder [23], the 
height of the barrier AW in Fig. 3 yields the line tension 
(j i = AW/2D, with D the short edge of the rectangle. For 
zb = 6.0 we obtain oi w 0.47 (in units of fcgT per particle 
length). As expected, the line tension increases rapidly 
with increasing z^, as manifested by the growing barriers 
of Fig. 1. 

In analogy to the liquid-gas transition, we can construct 
a binodal, by plotting the peak positions in P(N) as a 
function of Zb (Fig- 4). In the thermodynamic limit L — ► 
oo, the gas and liquid branches of the binodal meet at the 
critical point, while in finite systems they "sway" around 
it (finite-size rounding [24]). The horizontal line in Fig. 4 
marks the thermodynamic limit estimate of the critical 
fugacity Zb,ct, taken from Fig. 5. Also shown in Fig. 4 
are the "diameters" , defined as the average of the gas and 
liquid peak positions. In contrast to the binodal, finite size 
effects in the diameter arc much weaker, the reason being 
that the singularity in the latter is only logarithmic in 
two dimensional Ising systems [25]. From the intersection 
of the diameter with the line zg iCr in Fig. 4, we obtain 
p cr « 2.42 for the critical density of A particles. Due to 
symmetry, the overall particle density is twice this value. 

The critical fugacity zg cr was obtained from a scaling 
analysis of the finite-size susceptibility [22] 
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Fig. 2: Snapshot of the 2D Zwanzig model at coexistence in a 
rectangular L x D simulation box, with L = 30 and D = 15, 
at fugacity zb = 6.0. 




Fig. 3: W(N) as obtained in a rectangular 25 x 10 system at 
zb = 6.0. Note the flat region in between the peaks, and also 
the definition of the barrier AW. 



5.6 



5.5 



5.4 



5.3 



5.2 
















XL = 



(MY 



Fig. 4: Phase diagram (binodal) of the 2D Zwanzig model for 
system sizes L = 15, 20, . . . , 40 (from outer to inner). The left 
(right) branches mark the positions of the gas (liquid) peak, the 
horizontal line is the critical fugacity zb,ci ~ 5.294 obtained 
from Fig. 5(a). Also shown are the diameters (middle curves). 
The intersection of the diameter with the horizontal line yields 
the critical point (circle). 
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Fig. 5: Finite size scaling analysis of the susceptibility \l- (a) 
Positions of the peak maxima Zb,l versus 1/L. The line is a 
linear fit from which zb. cr follows, (b) The values xi,max at 
the maxima versus L. The curve is a fit to a powerlaw, see 
details in the text, from which the critical exponent 7 follows. 



with m = N— (N) . Following standard arguments [21] , xl 
versus zb exhibits a maximum, at position zb,l and with 
value xl, max- The peak positions scale with L according 
to zb,l — zs,cr oc L~ x l v , with v the critical exponent of 
the correlation length. Using the 2D Ising value v = 1, 
we observe excellent scaling, see Fig. 5(a), and by fitting 
we obtain zb.cx 

~ 5.294. Of course, due to symmetry, 
the critical fugacity of the A particles is also equal to this 
value. In addition, the peak maxima are expected to scale 
a s Xl, max oc LP 1 '■ ^ ', with 7 the critical exponent of the 
susceptibility. Indeed, the maxima scale accordingly, see 
Fig. 5(b), and by fitting we obtain 7 w 1.76, in good 
agreement with the 2D Ising value 72D.7 = 7/4. 

Discussion. — We have shown that the 2D Zwanzig 
model undergoes a liquid-gas transition at critical fugacity 
z cr w 5.294. As expected for systems with short-ranged in- 
teractions, the transition belongs to the universality class 
of the 2D Ising model. Note that the model studied here 
corresponds to that of hard rods on square lattices, in the 
limit where the rod length k — ► 00. As was shown in [6], 
the square lattice variant also exhibits a 2D Ising critical 
point, consistent with our findings. In these and other 
works [5-8], the resulting order above z cr is termed ne- 
matic, and the corresponding transition at z CI an isotropic- 
to-nematic transition. In contrast, our work indicates that 
the transition is just the liquid-gas transition. Therefore, 
the resulting order should be termed magnetic, since the 
liquid-gas transition is isomorphic to the formation of a 
spontaneous magnetization in the Ising model below its 
critical temperature. 



Further evidence in favor of a liquid-gas transition, and 
against isotropic-to-nematic, is the coexistence between 
ordered states, see Fig. 2. Note that phase coexistence in 
the 2D Zwanzig model is possible for all fugacities in the 
ordered region z > z ct of the phase diagram. The latter is 
analogous to the coexistence between domains of negative 
and positive magnetization in the Ising model below its 
critical temperature. In both cases, the line tension in- 
creases as one moves deeper into the ordered region. The 
corresponding interfaces are therefore order- order inter- 
faces, which (likely) do not exist between nematic phases 
of liquid crystals, where the particle orientations are con- 
tinuous [26,27]. Of course, liquid crystals may exhibit 
isotropic-nematic coexistence, at a first-order isotropic-to- 
nematic transition. In that case one has order-disorder 
interfaces, which survive only at the transition point. 

For the 2D Widom-Rowlinson (WR) model consist- 
ing of disks, the (single-species) critical density equals 
p CI ~ 0.78, while for the critical fugacity z CI sa 1.73 is 
obtained [16]. These values arc significantly below the 
Zwanzig values reported here, implying a rather weak de- 
pletion effect in the latter. This can be made plausible by 
considering the excluded volume per particle. In a fluid of 
line segments, unlike segments exclude a volume I , with 
I being the length. In a mixture of disks, one obtains irl 2 , 
with / being the disk diameter. Hence, the excluded vol- 
ume per particle is 7r times larger in fluids consisting of 
disks, and so we expect the transition at a density and 
fugacity reduced by roughly the same factor. Indeed, in- 
spection of the reported numerical estimates follow this 
prediction reasonably well (to within 3%). 

Our main conclusion is therefore that phase transitions 
observed in the 2D Zwanzig model and its lattice variants 
should not be compared to liquid crystal transitions, but 
instead to transitions observed in simple fluids. This not 
only includes liquid-gas transitions, but also the possibility 
of crystallization at high density. Of course, the model 
considered in the present work cannot crystallize, since 
the particles are infinitely thin, but the lattice variants 
considered elsewhere may [5-8] . Interestingly, these works 
indeed report evidence of a second transition occurring at 
high density; whether this transition can be interpreted in 
terms of (quasi) crystallization [28] could be an interesting 
topic. 

Finally, we discuss the expected trends when the 2D 
Zwanzig model is extended to a larger set of particle ori- 
entations. For hard rods on triangular lattices, the uni- 
versality class is that of the 2D 3-state Potts model [6]. 
We expect the same universality class for an off-lattice 
fluid of hard line segments, with three allowed orientations 
9 G {0, 7r/3, 27r/3} per segment (with 8 the angle between, 
say, the segment and the x-axes) . In the Potts model [29] , 
the nearest-neighbor spin interaction assumes two values: 
a low (energetically favorable) value when two neighboring 
spins are in the same state, and a higher value when they 
are not. For line segments, the analogue is the excluded 
volume ad = |sin^| between pairs of segments, with </> 
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the angle between the segments. The reader can verify 
that for the above set of three orientations, one either has 
ad = (when two segments are aligned) or = V3/2 
(when they are not). The excluded volume term thus ex- 
hibits the same symmetry as the Potts pair interaction, 
which naturally accounts for the observed critical behav- 
ior on triangular lattices [6]. However, for larger sets of 
allowed orientations, we expect the analogy to the Potts 
model to break down. For instance, allowing four orienta- 
tions 9 £ {0, 7r/4, 7r/2, 37t/4}, the excluded volume already 
assumes three values, in disagreement with Potts symme- 
try. Note that the analogy to the Potts model may also 
be useful to interpret results obtained in three dimensions. 
In this case, the Zwanzig model consists of hard rods of 
finite width, allowed to point in three mutually perpen- 
dicular directions. In terms of symmetry, this corresponds 
to a three-dimensional 3-state Potts model, which has a 
first-order phase transition [30,31]. The latter is indeed 
consistent with Zwanzig's result, where the existence of a 
first-order transition was also demonstrated [12]. In three 
dimensions, it should also be possible to realize 4-state 
Potts symmetry, by choosing the rod orientations perpen- 
dicular to the faces of a regular tetrahedron. For larger 
sets of orientations, the analogy to the Potts model is 
again expected to break down. We refer the interested 
reader to [32,33] where Zwanzig models with large sets of 
allowed particle orientations are discussed. 

In summary, we have presented simulation results of 
the 2D Zwanzig model. In agreement with lattice vari- 
ants of this model, we find a phase transition with critical 
point belonging to the universality class of the 2D Ising 
model. The novelty of the present work has been the in- 
terpretation of this transition in terms of the liquid-gas 
transition, as opposed to isotropic-to-nematic. This in- 
terpretation accounts naturally for the observed critical 
behavior, as well as for the observed phase coexistence 
in the ordered region of the phase diagram. In addition, 
the excluded volume interaction in fluids of hard line seg- 
ments with restricted orientations, was shown to resemble 
in some cases the symmetry of the g-state Potts model 
with q > 2, thereby elucidating a previous simulation re- 
sult [6]. 
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